home *** CD-ROM | disk | FTP | other *** search
/ Languguage OS 2 / Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO / language / dino / dino_bot.1 / examples / complex / lu.d < prev    next >
Encoding:
Text File  |  1991-03-10  |  3.8 KB  |  133 lines
  1. /* Copyright, 1990, Regents of the University of Colorado */
  2. /* This program computes the LU decomposition of a matrix, with partial
  3.  * pivoting.  The last row of A contains pivoting information - what
  4.  * row was swapped with at each iteration of the algorithm. */
  5.  
  6. #define P 6
  7. #define N 14
  8.  
  9. #include "dino.h"
  10.  
  11. environment node[P:id] {
  12.  
  13. #include <math.h>
  14. #include <stdio.h>
  15.  
  16. composite plu (a)
  17.  
  18. double distributed a[N+1][N] map WrapCol;
  19.  
  20.                                 /* ==> Declares the variable a to be mapped
  21.                                        by columns to the processors, but in
  22.                                        a wrap fashion, so that each processor
  23.                                        gets columns id, id+P, id+2P, ... */
  24.  
  25. {
  26.   int i, j, k;          /* Looping variables */
  27.   double piv;           /* Value of the pivot */
  28.   int pivrow;           /* Index of the pivot */
  29.   double temp;          /* Temporary */
  30.  
  31.   double distributed m[N+1] map all;    /* Holds the multipliers */
  32.  
  33.                                 /* ==> Illustrates the declaration of local
  34.                                        distributed data within a composite
  35.                                        procedure. */
  36.  
  37.   /* Loop through the passes of the algorithm */
  38.   for (i = 0; i < N; i++) {
  39.  
  40.     /* Check to see if we have the column for this iteration */
  41.     if (i % P == id) {
  42.       /* Select the pivot row */
  43.       piv = fabs (a[i][i]);
  44.       pivrow = i;
  45.       for (j = i + 1; j < N; j++)
  46.         if (fabs (a[j][i]) > piv) {
  47.           piv = fabs (a[j][i]);
  48.           pivrow = j;
  49.         }
  50.  
  51.       /* Swap in the pivot column */
  52.       temp = a[pivrow][i]; a[pivrow][i] = a[i][i]; a[i][i] = temp;
  53.  
  54.       /* Compute the multipliers */
  55.       for (j = i + 1; j < N; j++)
  56.         a[j][i] /= a[i][i];
  57.  
  58.       /* Put the pivot row into a, in preparation for sending */
  59.       a[N][i] = pivrow;
  60.  
  61.       /* Send out the data to the multiplier vectors on all processors */
  62.       m[<i+1,N>] = a[<i+1,N>][i];
  63.  
  64.                                 /* ==> Ranges may be used in assignment
  65.                                        statements to perform array
  66.                                        assignments. */
  67.  
  68.       m[<i+1,N>]# {node[] - node[id]} = m[<i+1,N>];
  69.  
  70.                                 /* ==> Explicit sends are indicated by putting
  71.                                        an environment set in braces after the
  72.                                        "#" sign. */
  73.  
  74.                                 /* ==> Sends to every node but itself through
  75.                                        the use of the "-" operator in the
  76.                                        environment expression. */
  77.     }
  78.  
  79.     else
  80.  
  81.       /* Receive the multiplier data */
  82.       m[<i+1,N>] = m[<i+1,N>]# {node[i%P]};
  83.  
  84.     /* Now, we're ready to do the elimination and pivoting at the same time */
  85.     pivrow = m[N];
  86.     for (j = id; j < N; j += P)
  87.       if (j > i) {
  88.         temp = a[i][j]; a[i][j] = a[pivrow][j]; a[pivrow][j] = temp;
  89.         for (k = i+1; k < N; k++)
  90.           a[k][j] -= m[k] * a[i][j];
  91.       }
  92.   }
  93. }
  94.  
  95. }
  96.  
  97. environment host {
  98.  
  99.   double a[N+1][N];                     /* Holds the matrix */
  100.  
  101.   void main () {
  102.  
  103.     int i, j;                           /* Looping variables */
  104.  
  105.     /* Set up the test matrix */
  106.     for (i = 0; i < N; i++) {
  107.       for (j = 0; j < N; j++)
  108.         a[i][j] = (i + 1) * (j + 1);
  109.       a[i][i]++;
  110.     }
  111.  
  112.     /* Output an initial message */
  113.     printf ("Starting LU decomposition...\n\n");
  114.  
  115.     /* Call the composite procedure */
  116.     plu (a[][])#;
  117.  
  118.     /* Printout the results */
  119.     printf ("Result is...\n");
  120.     for (i = 0; i < N; i++) {
  121.       for (j = 0; j < N; j++)
  122.         printf ("%6.1f", a[i][j]);
  123.       printf ("\n");
  124.     }
  125.  
  126.     printf ("\nWith pivots...\n");
  127.     for (j = 0; j < N; j++)
  128.       printf ("%.2f ", a[N][j]);
  129.     printf ("\n");
  130.  
  131.   }
  132. }
  133.